Model reconstruction of nonlinear dynamical systems driven by noise 
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An efficient technique is introduced for modei inference of complex nonlinear dynamical systems 
driven by noise. The technique does not require extensive global optimization, provides optimal 
compensation for noise-induced errors and is robust in a broad range of dynamical models. It is 
applied to clinically measured blood pressure signal for the simultaneous inference of the strength, 
directionality, and the noise intensities in the nonlinear interaction between the cardiac and respi- 
ratory oscillations. 
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Most natural and man-made systems are inherently 
noisy and nonlinear. This has led to the use of stochas- 
tic nonlinear dynamical models for observed phenomena 
across many scientific disciplines. Examples range from 
lasers p] and molecular motors , to epidemiology [3l , 
to coupled matter-radiation systems in astrophysics [^|. 
In this approach a complex system is characterized by 
projecting it onto a specific dynamical model with pa- 
rameters obtained from the measured time-series data. 
In a great number of important problems the model is 
not usually known exactly from "first principles" and one 
is faced with a rather broad range of possible parametric 
models to consider. Furthermore, important "hidden" 
features of a model such as coupling coefficients between 
the dynamical degrees of freedom can be very difficult 
to extract due to the intricate interplay between noise 
and nonlinearity. These obstacles render the inference of 
stochastic nonlinear dynamical models from experimen- 
tal time series a formidable task, with no efficient general 
methods currently available for its solution. 

Deterministic inference techniques p| consistently fail 
to yield accurate parameter estimates in the presence 
of noise. The problem becomes even more complicated 
when both measurement noise as well as intrinsic dynam- 
ical noise are present |10| . Various numerical schemes 
have been proposed recently to deal with different aspects 
of this inverse problem HH i,H[ll [H [13. A standard 
approach is based on optimization of a certain cost func- 
tion (a likelihood function) at the values of the model pa- 
rameters that best reconstruct the measurements. It can 
be further generalized using a Bayesian formulation of 
the problem [8|, Il0| . Existing techniques usually employ 
numerical Monte Carlo techniques for complex optimiza- 
tion or multidimensional integration [g tasks. Infer- 
ence results from noisy observations are shown to be very 
sensitive to the specific choice of the likelihood function 
. Similarly, the correct choice of this function is one of 
the central questions in the inference of continuous-time 
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noise-driven dynamical models considered here. 

In this Letter, we present an efficient technique of 
Bayesian inference of nonlinear noise-driven dynamical 
models from time-series data that avoids extensive nu- 
merical optimization. It also guarantees optimum com- 
pensation of noise-induced errors by invoking the like- 
lihood function in the form of a path integral over the 
random trajectories of the dynamical system. The ro- 
bustness of our technique in a wide range of model pa- 
rameters is verified using synthetic data from the stochas- 
tic Lorenz system. We also present the reconstruction of 
a nonlinear model of cardio-respiratory interaction from 
experimentally measured blood pressure signal. 

Let the trajectory x(f) of a certain A-dimensional dy- 
namical system be observed at sequential time instants 
t ,ti,... and a series y = {(tk,Yk)', k = : K} thus 
be obtained. These measurements can be related to the 
(unknown) system trajectory through some conditional 
probability distribution function (PDF) p Q [y|x(t)] giv- 
ing the probability of observing a time series y for a 
specific system trajectory x(t). If we assume that yfc 
has the same dimension as x(ifc) and the measurement 
errors y& — x(t^) are uncorrelated Gaussian random vari- 
ables with mean zero and variance e 2 , then we obtain 
Po (y\X) = Uto M[yi -x(t, ),ej], where X = {x(t*)}. 

Assume now that the underlying dynamical system is 
in fact nonlinear and stochastic, evolving according to 

x(i)=f(x) +£(*), (I) 

where £(t) is an additive vector noise process. We param- 
eterize this system in the following way. The nonlinear 
vector field f(x) is written in the form 

f(x)=U(x)c = f(x;c), (2) 

where U(x) is an N x M matrix of suitably chosen basis 
functions {U nm (x); n = 1 : N, m = 1 : M}, and c is 
an M-dimensional coefficient vector. An important fea- 
ture of (|2J) for our subsequent development is that, while 
possibly highly nonlinear in x, f (x; c) is strictly linear in 
c. Dynamical noise £(£) may also be parameterized. For 
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instance, if is stationary white and Gaussian 

(€(*)> =0, (^)f(0)=D^-t'), (3) 

then the N x N (symmetric) noise covariance matrix D 
fully parameterizes the noise. The vector elements {c m } 
and the matrix elements {D nn i} together constitute a set 
M. = {c, D} of unknown parameters to be inferred from 
the measurements y. 

In the Bayesian model inference, two distinct PDFs are 
ascribed to the set of unknown model parameters: the 
prior p pr (A4) and the posterior p pos t(M.\y) , respectively 
representing our state of knowledge about M before and 
after processing a block of data y. The two PDFs are 
related to each other via Bayes' theorem: 

i(y\M) Ppi (M) 



here we introduce the "velocity" y k and matrix $(x) 
Yfe = h^ 1 (y k+1 -yk), (*(x;c))„„' =df n (x.;c)/dx n >. 

With the use of J2J), substitution of the prior p pl (A4) 
and the likelihood £(y\M) into 0} yields the posterior 
P P ost(M\y) = const x exp[-S(M\y)}, where 



S(M\y) = S y (c,D) = i Py (D)-c T w y (D) + Vs y (D)c. 



Here, use was made of the definitions 

K-l 
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(I)) = hJ2 yl & 1 + K ln(det D), 



(7) 
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P V ost{M\y) 



J £(y\M) Ppr (M)dM' 
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Here £(y\A4), usually termed the likelihood, is the condi- 
tional PDF of the measurements y for a given choice M. 
of the dynamical model. In practice, Q can be applied 
iteratively using a sequence of data blocks y, y', etc. The 
posterior computed from block y serves as the prior for 
the next block y' , etc. For a sufficiently large number of 
observations, p poB t{M\y,y' , . . .) is sharply peaked at a 
certain most probable model M. = M* . 

We specify a prior distribution p pr (A4) that is Gaus- 
sian with respect to elements of c and uniform with re- 
spect to elements of D. Thus, p pl (M) = Af(c—c pI , E pr ), 
where the mean c pr and the covariance E pr respectively 
encapsulate our knowledge and associated uncertainty 
about the coefficient vector c. We now write the ex- 
pression for the likelihood in the form of a path integral 
over the random trajectories of the system: 



U^D- 1 y fc --v(y fc ) 



w y (D) = S pr 1 c pr + /i^ 

fc=0 
K-l 

S y (D) = S- 1 + / l ^U^D" 1 ^, 

k=0 

where = U(yfc) and the components of vector v(x) 
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>(*) = E ^r (x) . „, = 1: ^. (in 



i(y\M) = 



X(tf) 



dx r , 



The mean values of c and D in the posterior distribu- 
tion give the best estimates for the model parameters for 
a given block of data y of length K and provide a global 
minimum to Sy(c, D). We handle this optimization prob- 
lem in the following way. Assume for the moment that c 
is known in JJJ. Then the posterior distribution over D 



x(t,) 



Po(y\<%) J~m [ x (0] ^ x (^)i (5) has a mean D' t = y (c) that provides a minimum to 



where we choose t{ <C to < ^ U so that I does not 
depend on the particular initial and final states x(t;), 
x(tf). The form of the probability functional Tm over the 
system trajectory x(£) is determined by the properties of 
the dynamical noise 0,0]. 

In this Letter, we are focusing on the case of Gaussian 
white noise, as indicated in (JJJ, ©. We consider a uni- 
form sampling scheme t k — to + hk, h = (tx — to)/K 
and assume that for each trajectory component x n (t) 
the measurement error e is negligible compared with 
the fluctuations induced by the dynamical noise; that 



S y (c, D) with respect to D = D T . Its matrix elements 
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(c) = E \j k ~ U(y fe ) c 

fe=0 



Yk - U(y fe )c 



(10) 

Alternatively, assume next that D is known, and note 
from (JZJ) that in this case the posterior distribution over 
c is Gaussian. Its covariance is given by H y (D) and the 
mean c post minimizes S l y (c, D) with respect to c 



is, e 2 < h(D 2 



Uk=o^iyk — x (tk)} in ©. Using results from [Tj| for 
^m[x(£)], the logarithm of the likelihood © takes the 
following form for sufficiently large K (small time step 
h): 



Consequently, we use p Q (y\X) 
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H y - 1 (D)w y (D) 
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--bgf(y|M)=lndet6 + -E [tr*(y fc ;c) (6) 



k=Q 



(y fc -f(y & ;c)) J D 1 {y k - f(y fc ;c)) +N\n(27rh) 



We repeat this two-step optimization procedure itera- 
tively, starting from some prior values c pr and S pr . We 
do not need a prior for D at this stage, according to 
(|5|)- (|llfl . At convergence we obtain the "true" mean pos- 
terior values, c post — > c post and Dp OSt — > D post . The 

posterior covariance matrix X post = H~ 1 (D post ). 

To continue the inference process with a new block 
of data y' of length K' we update the prior mean 
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~pOSt ; -Dp r 



Dp OS t, and covariance, E pr 



J post ) 



and repeat the two-step optimization procedure. The 
modification is that updates now explicitly involve D pr 

K ~ K' 

:D pr 



D 



post 



K' + K 



K' + K 



&y,{C p0st ). (12) 



We obtain l|12|) from l|10|l where the data record y U 
y' of length K + K' is used instead of y and the sum 
over the first K data points (block y) is given by the 
D pr . Clearly, many non-overlapping, and not necessarily 
contiguous, data blocks of varying lengths may be used 
in this recursive model inference algorithm fla ]. 

The terms involving tr<t(y k ) in © originate from the 
prefactor in J-jvi [ x (^)] © an d do not vanish at the dy- 
namical system attractors {Q, unlike the terms in JHJl in- 
volving y k — f (ys,-; c). Therefore both types of terms are 
required to optimally balance the effect of noise effect in 
{yfc} an d provide the robust convergence. 

We now consider an example of the dynamics in l|T)l 
given by a noise-driven chaotic Lorenz system. It has 
dynamical variables, x = {xi, x 2 , x$}, the vector field 

f (x) = (#2 — %1, rxi - x 2 - xix 3 , xix 2 -bx 3 ), (13) 

and noise correlation matrix (£ n {t)£n' (£')) = d n 5 njn >. 
The parameters in Ijl3(l are a — 10, r = 28, b = |. 
We sample a system trajectory x(t) and produce a "data 
record" {y(t k )} to be fed directly into the algorithm. As 
an inferential framework, we introduce the following (bi- 
linear) model of stochastic dynamics for x(i): 



Xji ^ Q"ni Xi{t^j 
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Kj=l 



!,(*) &,(*), (14) 



n,i, j = 1,2,3, where M = {{a ni }, {b ni j}, {£>„„'}} is the 
vector of 18 unknown coefficients. The form of the 12 
basis functions U nm (x) is evident from 1|14|) . We were 
able to infer the accurate values of c, D for time step h 
varying from 0.01 to 10 -6 and noise intensities d n varying 
from to 10 2 . An example of convergence of the coeffi- 
cients is shown in Fig. ^ We found a step-wise decrease 
in variances that occurs on a time scale of the period of 
oscillations t osc » 0.6 (dashed line in the figure). The 
error of the inference is sensitive to the noise intensity, 
total time T and the time step h. For example, for the 
parameters of the curve 1 in the Fig. ^thc relative error 
was 0.015%. The ratio T/h has to be increased at least 
250 times to achieve error less then 1% when the noise 
intensity is increased 10 4 times (curve 2 in the figure). 

Finally, we apply our method to study the stochastic 
nonlinear dynamics of complex physiological system. To 
be specific we infer the strength, directionality and a de- 
gree of randomness of the cardiorespiratory interaction 
from the central venous blood pressure signal (record 
24 of the MGH/MF Waveform Database available at 
www.physionet.org). Such estimations provide valuable 
diagnostic information about the responses of the au- 
tonomous nervous system |l6l Il7|. However, it is inher- 
ently difficult to dissociate a specific response from the 




FIG. 1: Examples of convergence of the coefficient 0,31 cor- 
responding to parameter r of the Lorenz system for the 
total length of the time record T = 560 and the following sets 
of parameters: line 1 {d n } = {0.01,0.012,0.014} time step 
h = 0.002; line 2 {d n } = {100,120,140} and h = 0.00002. 
The insert shows dispersion (Ao|i)={a3x) - {«3i) 2 for the 
same sets of parameters. The vertical dashed line shows the 
time-scale of the step-wise decrease in the variance. 



rest of the cardiovascular interactions and the mechanical 
properties of the cardiovascular system in the intact or- 
ganism [lj| ■ Therefore a number of numerical techniques 
were introduced to address this problem using e.g. linear 
approximations [l9j . or semi-quantitative estimations of 
either the strength of some of the nonlinear terms [2(j 
or the directionality of coupling 0, . But the prob- 
lem remains wide open because of the complexity and 
nonlinearity of the cardiovascular interactions. Our al- 
gorithm provides an alternative effective approach to the 
solution of this problem. To demonstrate this we use a 
combination of low- and high-pass Butterworth filters to 
decompose the blood pressure signal into 2-dimensional 
time series {s(t k ) = (so(tfc), si(t k )), t k = kh, k = : 
K} representing observations of mechanical cardiac and 
respiratory degrees of freedom on a discrete time grid 
with step h = 0.002 sec. We now introduce an auxil- 
iary two-dimensional dynamical system whos trajectory 
x(t) = (xo(t), Xi(t)) is related to the observations {s(ifc)} 
as follows 

/. \ s n (t k + h)- s n (t k ) , , 

X n {tk) — 0,ln ; 1" a 2nSn\tk) + ^3 n , 

h 

where n = 0,1. The corresponding simplified model of 
the nonlinear interaction between the cardiac and respi- 
ratory limit cycles has the form (cf. with p^l 

x n = hn + b 2n Sn + h n x n + b in s 2 n + b 5n x1 + b 6n s n x n 
+ h n s n + b Sn s n x n + b 9n x n x n + b Wn x n + bu n x n x n -i 

+ bl 2n X n X 1 _ n + b 13n X n xl_ n + £ n (t), 71 = 0,1. (15) 

where £ n (t) is a Gaussian white noise with correlation 
matrix D nn > J3J). We emphasize that a number of impor- 
tant parameters of the decomposition of the original sig- 
nal (including the bandwidth, the order of the filters and 
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FIG. 2: (a) Time series of the cardiac oscillations S2(n) in ar- 
bitrary units (black line) obtained from central venous blood 
pressure with 4-th order Butterworth filter (low-and high- cut- 
off frequencies were /; = 0.8Hz and fh = 2.8Hz and the sam- 
ple rate was 90 Hz after resampling of the original signal). 
Inferred time series of the cardiac oscillator (green line), (b) 
Power spectrum of the cardiac oscillations obtained from the 
real data (black line). Power spectrum of the inferred oscilla- 
tions (green line), (c) Limit cycle of the cardiac oscillations 
(x2(n), ?/2 (n-) obtained from real data as described in the text 
(black line). Limit cycle of inferred oscillations (green line). 

the scaling parameters aui) have to be selected to mini- 
mize the cost {7J and provide the best fit to the measured 
time series {s(ifc)}. The parameters of the model (|15|) can 
now be inferred directly from the noninvasivcly measured 
time series of blood pressure. The comparison between 
the time series of the inferred and actual cardiac oscil- 



lations is shown in Fig. [21 Similar results are obtained 
for the respiratory oscillations. In particular, the param- 
eters of the nonlinear coupling and of the noise intensity 
of the cardiac oscillations are b\\2 — 3.9, 6122 = 0.62, 
6132 = -13.4, and D 22 = 4.75 ((£*,(<)) = D 22 ). Consis- 
tent with expectations, in all experiments the parameters 
of the nonlinear coupling are two orders of magnitude 
higher for the cardiac oscillations as compared to their 
values for the respiratory oscillations reflecting the fact 
that respiration strongly modulates cardiac oscillations, 
while the opposite effect of the cardiac oscillations on 
respiration is weak. Remarkably, our technique infers si- 
multaneously the strength, directionality of coupling and 
the noise intensities in the cardiorespiratory interaction 
directly from the non-invasively measured time series. 

In conclusion, we have derived an efficient technique 
for recursive inference of dynamical model parameters 
that does not require extensive numerical optimization 
and provides optimum compensation for the dynamical 
noise-induced errors. We verified the robustness of the 
technique in a very broad range of parameters of dy- 
namical models, using synthetic data from the chaotic 
noise-driven Lorenz system. Successful application of the 
technique to inference from real data of the nonlinear in- 
teraction between the cardiac and respiratory oscillations 
in the human cardiovascular system is particularly en- 
couraging, as it opens up a new avenue for the Bayesian 
inference of strongly nonlinear and noisy dynamical sys- 
tems with limited first-principles knowledge. Future ex- 
tensions will include more realistic observation schemes 
with "hidden" variables and finite measurement noise. 
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